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Abstract 

A compact difference scheme is described for treating the first-order 
system of PDE's which describe the equilibrium equations of an elastic body. 
An algebraic simplification enables the solution to be obtained by standard 
direct or iterative techniques. 
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INTRODUCTION 


The conditions for the static equilibrium of an elastic body are 
described by an elliptic system of nine partial differential equations for the 
displacements and stresses. This paper describes a finite difference scheme 
which can be solved by standard direct or iterative methods and yields a 
solution which approximates smooth displacements with second-order accuracy. 
Iterative techniques can be attractive as a means of solving three-dimensional 
problems because they minimize storage requirements and present an algorithmic 
structure well-suited to advanced computer architectures. For material 
problems, these features can be useful for solving layered composite materials 
as well as materials with nonlinear properties. 

As described here, a serious limitation of this method vis-a-vis finite 
element methods is that it is applicable only to bodies which can be 
subdivided into cube-like volume cells. However, a means of removing this 

restriction will be described in another paper. 

Part I describes an algebraic approximation to the equilibrium conditions 
in a cell as expressed by tractions and displacements on the surface of the 
cell. The condition that traction forces balance across cell faces leads to 
an algebraic condition for equilibrium between any neighboring cells expressed 
solely in terms of displacements. A finite-sum approximation to the work due 
to tractions leads to an energy estimate and to a variational description of 
the algebraic equilibrium equations. 

Part II illustrates this development for an isotropic material using a 
plane stress assumption to reduce the problem to two dimensions. Several 
simple iteration schemes are used to investigate the numerical convergence of 
the method when a singularity is present. 

The methods described in this paper are closely related to those 

described by the authors in the context of a simpler problem [1]. 
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PART I 

General Development 


1.1 The Equilibrium Problem 

In this section we describe the equilibrium equations for an elastic body 
and attempt to motivate the origin of the finite difference scheme which will 
be described in the following section. 

We consider a material body occupying a domain Q on whose boundary, T, 

- i5 3n ° Utward unit normal; u = ( Uj , u 2 , u 3 ) T denotes the displacement 
vector, x (Tj , t 2 , T^) the symmetric stress tensor, and e = ( e j> e 2 > G ) 
the symmetric strain tensor. 

In the absence of body forces, the equations of equilibrium are described 
by the following system of first-order partial differential equations: 

a) 3 r + 3 t + 3 t = o 

x : “I x 2 -2 x 3 -3 U ’ 

^ e = (grad u + grad T u) /2 in J2, 


:) t. = l c. . e. 
-i . L , ij -i 
J=1 J 


i = 1,2,3. 


Here (a) states the conditions for the equilibrium of forces, (b) defines the 
strain tensor m terms of the displacement, and (c) is the constitutive 
relationship between stress and strain (Hooke's Law). In (c), ( Cij ) involves 
21 parameters. By assumption 


e x = l e. . x. . > 0 

ij ij 


( 2 ) 
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with equality holding iff e = 0. 

On the boundary surface Y, we let p = x*n denote the traction due to 

the stress. Then boundary conditions associated with (1) are 

u = u^ on 

(3) 

0 r 

p = £ on ^2 9 

where and Y ^ form a disjoint partition of Y ; if T ^ = 4> th e 

solution of (1) will be determined to within a rigid body displacement. 

As a result of solving this boundary value problem, the tractions £ 

will be determined on T and the displacements u_ on Y^ . Thus the 
tractions on Y, say p(T), will be related to the displacements on u(T) 
on T, a fact we indicate symbolically by 

(4) p(r) = R q u(D. 

We call the boundary operator a transmission operator. 

Let II (ft) = {to} denote a partition of ft into volume cells. Standard 
integration arguments show that u must be continuous and the surface 
tractions must balance across cell faces. Clearly, a necessary and sufficient 
condition for equilibrium in ft is that any individual cell be in equilibrium 
with any neighboring cell. 

Consider a cell to whose volume is O(h^), where h is a representative 
length scale, and whose boundary surface Y consists of m faces. In 
equilibrium, the tractions p(Y) on Y will be related to the displacements 
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u(y) on Y by means of a transmission operator , i.e., p(y) = R u(y). 
Let < <J ) >y indicate a vector whose m components represent the average values 
of <f> on the faces of to. We may then relate the average tractions <p> on 
the faces of U) to the average displacements on the faces by 

(5) <£> r ' ■£ <HV 

where R^ is an mx m transmission matrix which is related to R . When 

03 

conditions for the balance of average surface tractions across cell faces are 
adjoined to (5), as well as boundary conditions for average tractions and 
displacements on T , we may expect that the resulting system of algebraic 
equations for average displacement values will yield an approximation to the 
equilibrium problem for (1) - (3) as h 0 . 

This approach can be made practical only if the transmission matrices 
R u can approximated without an a priori knowledge of the 

transmission operators related to the continuous problem. A general 

method for constructing R^ on arbitrary cells will be described in a 
separate paper by one of the authors. 

In this paper we describe a simple construction for R^ which results 
when cubical cells are employed. In this case, the equilibrium conditions (5) 
and the balance of traction conditions can be given a particularly simple form 
using finite difference notations and which we then call a compact finite 
difference scheme. Finite difference methods can then be used to obtain a 
finite sum energy estimate and also to characterize the average displacement 
values as the solution of a quadratic variational problem. As a result, the 
algebraic solution of the compact scheme can be obtained by either direct or 
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iterative methods using a priori facts about the structure of the algebraic 
system. This development is explained below. 


1.2 A Compact Difference Scheme 

In this section we describe a compact difference scheme (eq. (10)) which 
approximately describes the equilibrium of a cell due to displacements and 
tractions on the cell faces. We then develop an energy estimate and present a 
variational principle for the difference scheme. An appropriate method for 
solving the difference scheme is given in the following section. 

We suppose that ft can be partitioned into regular cubical cells whose 
faces are parallel to the coordinate axes; w(x) indicates a cell with 
centerpoint at x = (x^ x 2> x^. We denote the average value of a function 
<J) on a face whose centerpoint is K by ()>(£). If h^, = Ax^/2, i - 1,2,3, 
the volume of w is Aw = Sh^h^h^. 

Next, define central average and central divided difference operators 

and 6^ by 

Wj ♦(*) S 0^ + h 1? x 2 , x 3 ) + 4>(x 1 - hj, x 2 , x 3 ))/2 


( 6 ) 


<}>(x) i (<Kx, + h, , x 9 


V 


- ♦(* 1 " 


|))/2h, 


The operators lJ 2 , ^2 , ^3 anc * ^3 are ^ e ^ ne< * similarly. 
Also define 

grad h <|> = (6^ <|> , 6^ c|> 9 $) 


div h - " 6 1 U 1 + 6 2 U 2 + 6 3 V 


( 7 ) 
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Finally, we write 


(UT) = (Uj T,, U 2 Tj, ^ T 3 ) 


® - (£^ ’ £ 2 ’ £3 ^ * 


The restriction to cube-like cells simplifies the evaluation of surface 
tractions. Let 1 = 1,2,3, denote the centerpoints of the opposite faces 
of a cell to. The outward normals n(S*) satisfy n(^) = -n(S7) so that 
the average surface tractions j>(£*) are given by 


pUb = ± t.( 5 *), 


i = 1,2,3. 


In the following discussion the stress components arise as traction 
forces on the faces of cells. An important consequence is that then the 
conditions for the balance of average traction forces across cell faces simply 
reduce to the conditions that the jump in value of vanish across a face 

X£ = const. , i = 1,2,3. 

Corresponding to the equilibrium equations (1), we propose to consider 
the following compact scheme: in each cell w 


»> 5 1 II * 5 2 I 2 * I 3 " 0 


b) p. x. = y 

-1 1 t-« 


iJ -J 


i = 1,2,3 


c) e - (grad^ u + grad^ u) /2 


d) P. u - k 2 h 2 6 t. = y u - k 2 . h 2 . 6 . t., 
1 111-1 J - 


i * i 


i,j = 1,2,3, 
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Here indicates a positive diagonal matrix. Clearly, (a), (b), and (c) 

are consistent with (1); the significance of (d) will become clear from later 
developments . 

Figure 1 indicates, with reference to a rectangular cell, the variables 
associated with the sides of the cell by the scheme (10). 

If the strain tensor (e) is eliminated from these equations there thus 
results a system of 18 algebraic equations for the 18 components of _u_ and 
the 18 components of the tractions on the faces of a cell. 

We now indicate how (10) leads to a finite-sum energy estimate. Recall 
that the work W done on the body can be evaluated from (1) by the use of 
integration-by-parts and the use of Gauss' Theorem with the result 

TO TO T 

(11) 2W = £ / e . . t . . dw = / p u a y + / u p dy-/u div t do). 

i,j a 1J 1J t 1 r 2 a 

For finite differences summation-by-parts results from the identity 

(12) = ( U . <(0(6^ if/) + (jju if/)(6. ()>), i = 1 ,3 . 

using (6). Also, using (7), Gauss' Theorem holds in the form 

(13) £ div X ^ = ^ — Ay, 

03 cn h " f 

in which Ay is the area of a face of a cell on which n_ is the outward 


normal . 
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Using (12), 


(14) 


div h “ * “.I Hi * 6. t. + (e,T) h> 


where, from (10c), 


(15) 


(E ’ T) h =1 ii ( ^i Ii> 
1 = 1 


Next, using ( lOd) , 


I (ji. u )(6 T.) = X l 6. T. + l K. h 2 (6, T.) 

4 l—i — . l — l .lii—i 

1 l l 


where \ is some constant. Thus, if we define 


(16) 


t e »' c k = 


(e,T) + Y k? h?(6. t.) 2 , 

n ? 1 1 i — i 


then (14) assumes the form 


(17) 


div, u t = [e,t] + X 1 div^ 


v h - 


Recalling (2) and (10b), we see that (e ,-r ) h is a positive definite 

function of the strain tensor; hence, from (16), so also is [ e , T ] . Also, 

T 

recalling (9), div^ u t is seen to represent the work per unit volume done 
by traction on the faces of a cell 


to. 
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Summing (17) over cells in ft 

(18) 1 div u T xAu) = £ t e > T l h Aw + I A. T div^xAto, 

weft weft weft 

so that for a solution of the compact scheme (10) 


(19) 

l 

w eft 

div u T tAw = £ [ e , t ] h Ato 

weft 

V 

o 

• 



Also, 

using Gauss' 

Theorem in the form (13), 

(18) can 

be written 

(compare 

(ID) 






(20) 

l l*.U h 

weft 

Aw = ^ £ T u® Ay + ^ u T p^ 
r i 

Ay - l 
weft 

T 

X div^xAw, 


where 

use has been 

made of (9). In this 

equation, 

u 0, p 0 are to be 


interpreted as the average values on cell faces on T of the data given by 
(3). 

Consider a problem in which u_ is prescribed everywhere on T. The 

preceeding discussion can be used to verify that the average values of 
displacements on cell faces which solve the compact scheme (10) also solve the 
variational problem 

(21) min l [e,x] Aw = £ £ T u° Ay 

u weft T 

for u satisfying the boundary conditions on T. The Euler conditions for 
this problem simply express the balance of traction forces, expressed in terms 
by the use of (10), across cell faces. 


of u 
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As mentioned earlier we will describe elsewhere how these ideas may be 
adopted to treat cells having general shapes. For this reason we shall not 
present here the details of the convergence argument as it applies to (10) 
except to cite the result: the solution u_ h of (10) converges to the 

solution u_ of (1) with accuracy 0 ( h 2 ) in an J^-norm while p h ( u h ) 

converges to _p(u) with an accuracy 0(h). These remarks apply, of course, 
only to sufficiently smooth solutions of (l). 

An important feature of the compact scheme (10) is that it only employs 
values of the displacements and tractions which arise as average values on 

cell sides. This is in contrast to many finite element methods which employ 
edge and vertex values of u. 


1.3 Solution Method 

The compact scheme (10) may be solved by direct algebraic techniques such 
as Gauss elimination, considering the displacements and tractions as unknown 
variables. A preferable approach, which we shall now describe, is to 

eliminate the traction variables so as to obtain an algebraic system involving 
only the displacement variables. We first indicate, in general form, the 
steps which lead to this elimination. The specific result upon which 
numerical calculations can be performed is given by equation (32). 

Let Y.(co), i = 1,2, •••,6 indicate a face of a cell o> and write 


Mh (- V,u(y 2 ),..., u(y 6 )] t , 

t£] y = [P<Y 1 ) , p( Y 2 ) , • • • , p(Y 6 )] T , 


( 22 ) 
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where, for brevity, reference to (i) has been omitted. 

Eliminating the strains from (10) and recalling (9) we may solve (10) for 
the tractions [p]^ terms °f the displacements [u]^ to obtain 


(23) 


[p] = R [u] , 

— J y ur— J y 


where K is a block 6x6 transmission matrix associated with w which we 
0) 


write in terms of its rows as 


(24) 


R 


r h (Y,) 

— CO 1 


tl, N 


h, > 

I»<V 


(A direct evaluation using (10) shows that this matrix is symmetric.) Thus 
(23) states 


(25) p(Y.(a))) = r^(y i (o)))[u] Y((o) , i = 


rp 

As noted earlier, div, u t represents the work per unit volume due to 

h — 

the tractions arising from the displacements on the faces of u). Let 


Ay = diag(AY 1 ,AY 2 


•••,ay 6 ) . 
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In terms of the notations just described and recalling (19) we then have 


(26) 


0 < (div h u T t)Auj = [u] y T Ay [p] = [u]* Ay R^u] 


or y 


This shows that is positive semidef inite since the equality in (26) 

holds, according to (19), iff (pe) = 0 and this is seen to hold 

iff u(yJ = constant, i = 1,2, ••*,6. 

• * T 

Since [e,r]^ = div^ u t, the variational principle (21) also applies in 
the form 


(27) 


min y [u] T , . Ay R h [u] . . 

u WE Q " Y(u,) “ “ Y(«) 


The conditions for a minimum resulting from this problem are, as remarked 

a A 

earlier, simply the balance of average traction forces across any face y 
common to any neighboring cells which, using (25), may be written as 


(28) 


^ <Y) 1 H't( u ) * V (Y) tHl Y <u') ■ °- 


If y lies on , then also 


(29a) 


u(y) = u^(y) 


while if y lies on r , then 


^ <t) 1 hI y( „) - E°<r) 


(29b) 
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Earlier remarks show that this system has a unique solution for u(y) 
unless r = 

We now give the explicit result of applying (28) on cell faces. Using 

(10c) 


(30) 


e E ( £i’£ 2 ’£3 ) 


< 5 ^ ^ fi 2 U l + 6 1 U 2^ 2 ^ 3 U 1 + 5 1 U 3^ 2 


(6 2 U 1 + 6 1 U 2 )/2 6 2 U 2 (6 3 U 2 + V 3 )/2 


^ 3 U 1 + 5 i u 3^ 2 ( <5 3 u 2 + 5 2 U 3^ 2 5 3 U 3 


Let 


P i 


2 i 2 

K. h. , 

l l 


i = 1,2,3, 


(31) 


A = 


P 1 P 2 * P 2 P 3 * P 3 P 1 ’ 


and let U). . indicate a cell whose centerpoint is x = ( iAx, , jAx 9 , kAx~ ) , 

1 , j , k l z J 

Then across the faces of 0). . incident with oj. - . , , w . . . , , and 

i,j,k 1+1 , j , K 1,J+I,K 


w. . , , 

i, j,k+l 


we find 


(32a) 


1 l C H U * h l p 2 (“l ‘ "3^ * h l p 3 t U l ' “2^ 


OJ. . , 


1 l C H U - h l p 2^ p l ■ "3^ - h l p 3 ( u l - u 2^ 


03 . 


i+1 * 
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(32b) 


4 l C 2 1 -l * h 2 p 3 ^2 " “lb * h 2 M u 2 " "ik 


jtu. . 

ijJ»k 


A ] C 2J l h ‘ h 2 P 3^2 ' " h 2 P 1^2 ' 


w. 


i, j + 1 ,k 


(32c) 


A l C 3£ U + h 3 P 2^3 ~ v lh + h 3 p i^ 3 “ U 2^ 


to. . 

i, J»k 


A l °3Z -Z ~ h 3 p 2^3 " P l^ ” h 3 P l^ y 3 " W 2^ 


(A) . 

ij J>k+1 


In these equations we have left unspecified nine parameters which arise 
in the matrices K. = h. . These may, as indicated in the next section, 
be chosen for convenience (cf. [1]). 

We call (32) the stress-eliminated form of the compact scheme (10). 

In summary, then, the system of equilibrium equations (28), (29), in 
which is the transmission matrix for the compact scheme (10), arises as 
the Euler equations for a related positive definite quadratic variational 
problem. The system (32) is thus seen to be solvable by direct elimination 
methods without pivoting; if iterative methods are considered Gauss-Seidel and 
SOR methods are applicable. 
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Part II 

Example of an Isotropic Material 
II. 1 Equilibrium Equations 

We consider an isotropic material characterized by Young's modulus E 
and Poisson's ratio v. We use a plane stress assumption to formulate the 
problem in two dimensions, the ^ plane say. This involves setting 

= 0 and assuming that the stress components T H ’ T 22 ’ X 21 are 
independent of (Timoshenko and Goodier [2])* In two dimensions, the 
stress-strain relationship (lc) may be written in component form as follows: 


T 11 ■ "ll * " £ 22 


T 22 ° " S U + C £ 22 


T 12 “ T 21 0 £ 12 


or, in terms of the displacement _u, 


T 11 ' C U 1 * " 3 x., u 2 


T 22 ’ " \ U 1 * 5 3 Xj U 2 


T 12 ‘ T 21 ■I°( 3 x. “l * U 2 )- 


( 33 ) 
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The parameters 5, rj and a are given in terms of Young's modulus and 
Poisson's ratio as follows: 

„ E Ev E 

C = 5- , n = , a = . 

(1 - v ) (1 - v ) (1 + v) 

The quantity ^ a is known as the shear modulus. 


II. 2 Method of Solution 

In this section we write down the compact scheme for the two- 
dimensional case when square cells are employed and obtain the transmission 
matrix which relates the tractions and displacements in each cell. The 
properties of the resulting system are then discussed in the context of its 
iterat ive solution. 

Analogous with (10), and upon elimination of the strains, we have the 
following compact scheme for the components of displacement and stress: 


(34a) 



> 


“x, T 11 = ? U 1 * n % U 2 


\ T 22 = " S x, U 1 * 5 % V 


\ hi ’ % T 12 ‘ J "l * 5 Xj u 2 ) 


(34b) 
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K " M x \ = &1 ^ \ T 12 ‘ ° 2 ^ 6 x, T H 
x 2 1 2 




(34c) 


fy - y )u, = b li 6 t. 
^ x 2 X 1 2 *> 


2 ' 2 f - d 2 h 2 6 T 

22 Xj 21 


where k* - diag(a 2 , b 2 ), k \ - diag(c 2 , d 2 ) in the rotation of (10). Using 
<34 a) and patting 6, - a 2 * e 2 , 6 2 = b 2 * d 2 , (34c) become 


fy ~ V )u. = 0 . h 6 t 

^ x 2 x i 1 1 2 


12 


(W X 2 ' y x^ U 2 = 8 2 h \ T 22‘ 

Using (9) to eliminate the strains from (31) we may solve for the 
tractions [p] y in terms of the displacements M y in each cell 0) 

obtain 


(35) 


lElf ■ Va’f ’ 


where 

^ T 

Ip] = (p(?j), £ ( V’ £ (5 2^ ’ 

[u] y = (u(£;^), u(C 2 ) 5 n( 5 1 ) > H^2^ ’ 

h 

and in which the transmission matrix R y is given by 
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( 36 ) 


R h = 
(0 


Here 



1 + 0c n 

-d-0c 12 ) 

1 - 0C 11 

-(I * Sc 12 

0- 1 

-d -0c 21 ) 

I + 0C 22 

"(I + 0c 21 ) 

X - Oc 22 

2h 

1 - 0c n 

-(I + 0c 12 ) 

1 + 0c n 

-<I - 0c 12 


"(I + 0c 21 ) 

1 - ° C 22 

-(X - 9c 21 ) 

I - o= 22 


'11 


C 0 

0 V 2 a 


'22 


V 2 a 0 

0 c 


T 

°21 C 12 


0 

V 2 cr 


6 1 0 
0 6 


2 J 


h 

The matrix is obviously symmetric and is easily shown to be 

positive semi-definite for positive and 0^ Balancing the tractions 

vertical and horizontal faces common to neigboring cells, we obtain, 
with reference to Figure 2, 

(37a) 2(1 + 0 c n )u(P 0 ) + (I - 0 0 ^( 0 ^) + u(P 2 )) 

= (I - © c 12 )(u( Ql ) + u(Q 4 » + (I + 0 c 12 )(u(Q 0 ) + u(Q 3 )), 

(37b) 2(1 + 0 c 22 )u(Q 0 ) + (I - 0 c 22 )(u( Ql ) + u(Q )) 


(1-0 c 21 )(u( Pl ) + u(P 4 )) + (I + 0 c 21 )(u(P Q ) + u(P )), 
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These equations correspond to the stress-eliminated form (32). If 
r 2 = <(), i.e., the values of the displacements are prescribed on the boundary, 

then the coefficient matrix of the above problem is symmetric and positive 
definite. These properties have an important consequence so far as the SOR 

method is concerned, since for such a system the convergence of SOR is 
guaranteed for any value of the relaxation factor a in the range, 0 < a < 2, 
(see Varga [3], Young [4]). Similar conclusions can be expected when mixed 
boundary conditions are considered, 

A natural iterative method for solving (37) is line SOR. The method 
involves solving block tridgiagonal systems firstly along all horizontal lines 
and then along vertical lines. In the next section we will compare this and 
point relaxation methods to a simple problem. 


II. 3 Numerical Example 

Consider a square domain on whose vertical edges the displacements are 
fixed and which experiences a uniform load along its top horizontal surface. 
This situation is illustrated in Figure 3. 

The boundary conditions are 


u 

= (o,o) T 

on 

X = 

I 2 

H 

/-"V 

0 

0 

V-/ 

It 

on 

y = 

1 H 
ro 

= (0,-l) T 

on 

y = 


Our experiments were performed using the plane stress approximation with 
values of Young's modulus and Poisson's ratio given by E = 10 7 , v = 0.3. 
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The equilibrium displacement uj is anti-symmetric and U 2 is symmetric 
about the line x = V 2 . However, in the computations that were performed, no 
advantage was taken of this symmetry. 

For the choice of parameters = ^2 = an important simplifi- 

cation occurs in (34), viz. 

I - 0A = d iag( 1 - 2e/a, 0) , 

A 

I - 0A = d iag( 0 , 1 - 2 c/a). 

It is with the above values of 0^ and 0^ that the results in this section 
were obtained. 

In our experiments we used the point Gauss-Se ide 1 , point SOR. and line 
SOR methods. The parameters for the SOR methods were chosen to be the optimum 
ones for Laplace s equation. Initially, the value of u was taken to be zero 
at all the grid points. Keeping in mind that the method can be expected to 
yield only second-order accuracy, the iterations were terminated when the 
norm °f the residuals was less than 10”^. 

In Table I we show the dependence of the number of iterations required 
to attain the convergence criterion on the mesh size for various iterative 


schemes . 
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Table I. Dependence of Number of Iterations on Mesh Size 


h 

pt. Gauss-Seidel 

pt. SOR 

line SOR 

1/2 

42 

34 

— 

1/4 

100 

54 

30 

1/8 

325 

97 

51 

1/16 

1240 

191 

99 


These results indicate that the rates of convergence 
Seidel and SOR methods on this problem are 0(h 2 ) and 0(h) 
Table II contains the values of the £ 2 _norms the solutions 
eliminated equations (37) for different grid sizes. 


of the Gauss- 
respect ively . 
to the stress- 


Table II. 


—Norms of Numerical Solution 


h 

llljlj 

11 u 2 “ 2 

1/2 

0.0221 

0.0979 

1/4 

0.0086 

0.0710 

1/8 

0.0055 

0.0603 

1/16 

0.0048 

0.0568 

1/32 

0.0046 

0.0556 


rs ■ 

We note that the convergence of llu^D^ is 0(h ) while that of 
II u 2 II is 0(h 3 ^ 2 ) as h + 0. This degradation in behavior is due to the 
singularities which are located at the top corners. 
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The values of the tractions are calculated by substituting the values 
of u_ in (35). The components of displacement and traction display the 
relevant symmetric or anti-symmetric properties to four decimal places. There 

are two integral checks that can be carried out to verify the computations, 
namely 


1 1 

“ / T ii(0,y)dy + / x (l,y)dy = 0, 

0 11 0 11 

(35) 

/ T (1 ,y)dy = - I . 

0 Z 

These conditions express the integral condition 

/ x * n dy = 0 

r ~ 

using the assumed symmetry of t J 2 about x - I . These integrals were 

computed numerically using the midpoint rule. The first integral check held 
exactly while the second one was found to be correct to the number of decimal 
places specified in the displacement calculation. 

Various plots are given of the solution. Figures 4, 5, 6 are contour 

plots of the stress componets T n , , t J 2 respectively. These were 

obtained using 64 cells in each direction. Figures 7a and 7b show the 
principal stress vectors within each computational cell. The principal stress 

directions are defined to be those vectors x * 0 which satisfy the 

eigenproblem 


(t - Al)x = 0. 
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Concluding Remarks 

The treatment of equilibrium with a volume force _f requires no essential 

modifications to the method. In this case (la) has the form £ 3^^ - _f; 

i i 

correspondingly (10a) is modified to 


(10a' ) 



We leave the details of the consequent developments to the reader. 
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Figure 5. Plot of the Stress component t 
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